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We present a new scheme for the construction of highly localized lattice Wannier functions. The 
approach is based on a heuristic criterion for localization and takes the symmetry constraints into 
account from the start. We compare the local modes thus obtained with those generated by other 
schemes and find that they also provide a better description of the relevant vibrational subspace. 



I. INTRODUCTION 

Although the translational symmetry of a crystalline 
solid imposes a delocalized basis of Hamiltonian eigen- 
states (Bloch's functions), it is sometimes advantageous 
to consider a transformation to a new set of basis func- 
tions with a local character. Beyond the mathematical 
equivalence (both sets span the same space of states), 
a local viewpoint is better suited for the analysis of 
concepts such as bonding which are eminently local in 
character. Recent work on electronic Wannier functions 
has shown the usefulness of a local representation iruthe 
chemical characterization of a given band subspacep i 
the analysis of bonding topology irua disordered system, 
and in more formal developments. u 

The lattice dynamical problem is formally very similar 
to the electronic one: a set of Bloch eigenstates (normal 
modes) represents the collective vibrations of the atoms 
in the crystal. A basis change to a set of local displace- 
ment patterns (lattice Wannier functions or local modes) 
can in principle be achieved. So far, the main application 
for these local modes has been in the field of structural 
phase transitions. Typically, the behavior of a given dis- 
persion branch or set of branches determines the essential 
instabilities of the system, and the associated degrees of 
freedom enter into the construction of an effective Hamil- 
tonian which reproduces the relevant physics. Through 
the use of a localized basis set, the number of coupling 
terms in the effective Hamiltonian can be relatively small, 
easing the statistical mechanical treatment and the inter- 
pretation of the results, fn particular, the anharmonic 
terms in the effective Hamiltonian can be kept local (on- 
site) , in contrast with what happens in a reciprocal-space 
description. 

This local mode approach has been used extensively in 
the past to gain an understanding of the behavior of com- 
plex systems, but until recently the local variables were 
treated as dummy degrees of freedom in a semi-empirical 
model, with their interactions fitted to reproduce the ob- 
served phenomena, fn the last few years, a new approach, 
in which the effective Hamiltonian is parametrized on the 
basis of first principles calculations, has had great success 
in studies, of the phase transition sequences in perovskite 
oxides, an Central to the parametrization process is the 
explicit construction of lattice Wannier functions, and 



two schemes have beea proposed to carry it out. Zhong, 
Vanderbilt, and Rabea (ZVR) used the structure of the 
zone-center soft mode in perovskite BaTiOs to construct 
symmetry-adapted highly localized local modes. Sub- 
sequently, Rabe and WaghmareCI (RW) generalized this 
approach to reproduce the normal modes at several (typ- 
ically, high symmetry) points of the Brillouin Zone. 

While both the ZVR and RW approaches have been 
broadly successful in the specific problems for which 
they were conceived, in this paper we will argue that 
they are not completely satisfying in some respects. We 
will present a new procedure to generate lattice Wannier 
functions, an approach which makes use of the available 
symmetry information, produces local modes with a high 
degree of localization, enables a systematic improvement 
of their quality, and is straightforward to implement. 



II. METHOD 

We are interested in describing a relevant subspace 1Z 
of the full 3iVp-dimensional configuration space of a crys- 
tal with p atoms per unit cell. Typically, we can choose 
7Z as a complete band of dispersion branches (complete in 
the sense that it is invariant under the action of the space 
group of the crystals) . Associated to a branch j is a set 
of normal modes (3A^p-dimensional vectors) which 
are eigenvectors j-af the Fourier transform of the force- 
constant matrix.E3 (The displacement in the a cartesian 
direction of the atom k in cell 1 is given explicitly by 
UjQ.,K,a).) The normal modes transform according to 
irreducible representations of the little groups G k . These 
representations, considered over the whole BZ, determine 
the band symmetry. The relevant subspace 1Z is spanned 
by all the {u^} in the band, but it is clear that any trans- 
formation 



(i) 



will lead to a new basis of extended states which we will 
call Bloch modes. Here n is the band dimension, the 
number of dispersion branches in the band. 

Having thus specified the relevant subspace by means 
of Fourier space variables {u^}, the problem we tackle is 
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the construction of a new basis {w™} which is local, as 
opposed to extended, in character. Mathematically, the 
k label should be exchanged by a local label n associated 
to the different unit cells in the crystal. Translational 
symmetry takes the form: 

wj(l,K,a) =wf +t (l + t,K,a) , (2) 

which is trivially satisfied by the standard Wannier func- 
tion formw 

w l = 77 / cxp(-ikn)u k dk (3) 
" Jbz 

in which is the volume of the BZ. A high degree of lo- 
calization means that the displacement u>"(l, k, a) should 
be very small or zero when 1 is a few lattice constants 
away from n. The arbitrariness implicit in their defini- 
tion (Eq. 0) means that the Wannier functions are non- 
unique, and a relatively large latitude then exists to tune 
their properties. In particular, the degree of localization 
has traditionally been the focus-,of great interest, and re- 
cently, Marzari and Vanderbilttl have succeeded in opti- 
mizing the matrices appearing in Eq. to construct very 
localized electronic Wannier functions starting from the 
Bloch states. A restriction to unitary matrices resulted 
in an orthonormal basis of Wannier functions, and the 
optimization process led to symmetric-looking functions, 
even though no symmetry conditions were explicitly im- 
posed 113 In principle, such an approach should work for 
the vibrational problem, too. However, we prefer to take 
an alternate route which takes advantage of the knowl- 
edge of the band symmetry. 

A. Symmetry requirements 

As studied extensively in the literature, El one should 
supplement the translational constraints of Eq. |^ with 
another set of conditions which represent the transfor- 
mational properties of the wf under the effect of the 
point symmetry of the crystal. These are most easily 
discussed by introducing a symmetry-based definition of 
the center of a mode. Consider a Wyckoff set with rep- 
resentative site rj_and the set G T of operations in G that 
leave r invariant o Given an irreducible representation r 
of G r with dimension d T , any d T displacement patterns 
transforming with r under the action of G r are said to 
be centered in r. It is then notationally more convenient 
to use a double index to label these patterns: w r ^ s where 
s ranges from 1 to d T . The action of the elements of the 
space group G on this set generates images at the rest of 
the positions in the Wyckoff set, i.e., d r d T patterns u£. s 
per cell, where i ranges from 1 to d r (the multiplicity of 
the Wyckoff set). This set of lattice functions is repre- 
sented by the pair (r, r) and define a._Eepresentation of 
G which is called band represe.ntatian.E3 

A necessary condition for the description of a relevant 
band subspace by means of these symmetry-adapted local 



modes is the equivalence of the band symmetry of 1Z and 
the band representation (r, t) (in particular this implies 
n = d r d T ). More details about the choice of the correct 
(r, t) for a given 1Z are presented in the Appendix, where 
we also discuss the transformation properties of the cor- 
responding {Uj}- Incidentally, since Eq. || establishes a 
correspondence between lattice Wannier functions and 
Bloch modes, in what follows the latter can be also la- 
beled by the site and representation indexes: u k s 

B. Practical criterion for localization 

A straightforward scheme to obtain lattice Wannier 
functions can be based on a direct use of Eq. [| per- 
forming the BZ sum by means, of any of the standard 
"special k-points" methods.E3il3 The quality of the sub- 
space description can thus be systematically improved by 
simply using denser k-point sets. This approach can in- 
corporate information about the normal modes through- 
out the whole Brillouin zone, as opposed to at just one 
point (as in the ZVR method) , or at a very special set of 
high-symmetry k-points (as in the RW scheme). 

As stated in the Introduction, it is highly desirable 
that the local mode basis functions be as localized as 
possible, in order to permit the consideration of only a 
few coupling terms in the effective Hamiltonian. From 
the point of view of real applications, a basis of Wannier 
functions which are not localized is not efficient, even if 
it spans 1Z perfectly. The form of Eq. || suggests a very 
simple heuristic criterion to achieve a high degree of lo- 
calization for the lattice Wannier functions: choose the 
M k matrices in such a way that the u k Bloch vectors at 
different k add their contributions coherently at the cen- 
ter of the Wannier function. Interference effects can then 
be counted on to automatically dampen the amplitude of 
the displacements at sites away from the center. 

Both the symmetry requirements and the localization 
condition can be formulated in the following way. As- 
sume a (r, r) pair has been determined on the basis of 
band symmetry, and that we focus on the construction 
of local modes at cell n=0. Consider a set of d T (3Np- 
dimensional) orthonormal vectors {x r , s } which are cen- 
tered in r, transform with irreo-j and involve atoms in 
an orbit as close to r as possibleO The localization crite- 
rion is implemented by requiring that u k s be orthogonal 
to x Ti t if s ^ t, and "parallel'llmeaning that their scalar 
product is positive) if s = tS3 It can be seen that this 
condition fixes the form of the M k matrices, up to overall 
normalization factors. In general, the M k will not be uni- 
tary, with the result that two lattice Wannier functions 
at different cells w" s and s , will not beorthogonal if 
the pairs (r, s) and (r', s') are not equahEl In the next 
section we will provide a simple worked example of the 
new construction scheme and will compare its results to 
those of other methods. 
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III. EXAMPLES AND DISCUSSION 

In order to illustrate the scheme presented in the pre- 
vious section, we will employ a two-dimensional model 
crystal with two different atoms which occupy the la 
(white) and 16 (black) Wyckoff positions of the plane 
group p4mm (See Fig. [I] a)). A simple harmonic model 
for the force constants (with the couplings among the 
white atoms considered up to fourth nearest neighbors 
and the rest to first nearest neighbors, which corre- 
sponds to 6 independent parameters) gives the dispersion 
branches of panel b) in the figure. We will focus our at- 
tention on the two optical branches, which form a single 
band since they are essentially degenerate at the T and 
M points. These optical branches transform according to 
the decompositions 

T (4mm) : E 

X (2mm) : B x +B 2 (4) 
M (4mm) : E , 

in irreducible representations of the little co-groups at the 
high symmetry points. A simple application of the pro- 
cedure spelled out in the Appendix shows that the band 
representation compatible with the above band symme- 
try is that represented by the pair (o, E), in which E is 
a two-dimensional irreducible representation which turns 
out to be the vector representation of the point symme- 
try group at o. The set of {x} vectors is then trivial to 
construct: as the "o" Wyckoff position is occupied, it is 
just enough to make x .i and x ,2 unit vectors attached 
to the central atom and pointing in the x and y carte- 
sian directions, respectively. For this crystal structure, 
the simplest non-trivial set of special k-points is given 
by {(1/8, 1/8); (1/8,3/8); (3/8,3/8)}. The explicit appli- 
cation of the localization criterion proceeds as follows. 
At each k-point in the set the normal modes are com- 
puted and the M k matrices constructed. For example, 
at (1/8,3/8), the normal modes are 

u\= (...; 0.23, 0.93;...) 

i£= (...; 0.84, -0.19;...) , [0) 

where the ". . ." refer to displacements on atoms other 
than the one at the center. The "coherent addition at 
the center" condition then becomes: 

0.23 Mh +0.84 Mia > 

0.23 M 2 i + 0.84 M 22 = 

0.93 M n -0.19 M12 = [ ) 

0.93 M 2 i -0.19M 22 > , 

and is satisfied by 

/ 0.200 0.980 \ 

^0.964 -0.264 J ' { ' 

uniquely defined but for row-specific arbitrary factors. 
Since M is not unitary, the two optical Bloch vectors u 0jS 



at this k-point will not be orthogonal (although they can 
of course still be chosen to be normalized) . 

Once this procedure has been performed at every k- 
point in the set, the integral (sum) in Eq. [|can be carried 
out to give the components of the lattice Wannier func- 
tions. Since the {u k } determined by the localization cri- 
terion also satisfy the symmetry compatibility relations 
(see Appendix), the Wannier functions are symmetry- 
adapted. In Fig. U c) we show the displacements associ- 
ated to the local mode iu ,i, which transform as the first 
component of the vector representation (E) of 4mm. The 
degree of localization of these lattice Wannier functions 
can be gauged by computing the contribution to the to- 
tal norm from a given shell around the center atom, as 
presented on Table [|. Less than one per cent of the norm 
is outside the fourth shell (which corresponds roughly to 
the second- neighbor unit cells) . If the integral in Eq. || is 
computed using a denser special-point set, the degree of 
localization is maintained, as can be seen by comparing 
the columns labeled "this work (3 k)" and "this work (10 
k)" on Table B This means that the quality of the local 
modes can be systematically improved while retaining a 
high degree of localization. 

It is enlightening to. compare this scheme to that of 
Rabe and Waghmare.0 In the latter the analysis of the 
symmetry compatibility relations proceeds in the same 
way, and once the right (r, r) set has been identified, 
a series of orthonormal x sets is constructed at succes- 
sive shells centered on r. The extent of the outermost 
shell fixes the localization of the Wannier functions by 
construction, and the actual atomic displacements are 
determined by fitting to the normal modes computed 
at a few high-symmetry points of the Brillouin Zone. 
In essence, the normal-mode information determines the 
weight assigned to each symmetry-adapted shell, so there 
is a tradeoff between the extent of the lattice Wannier 
functions and the amount of information from the real 
dispersion relations that can be used in the construction 
procedure. For example, in the PbTiC>3 work, Rabe and 
Waghmare found that adding information about the nor- 
mal modes at the X point resulted in a less localized local 
mode than if only the (r, M, R) set was usedB In contrast, 
our scheme can deal with the extra k-point without loss 
of localization: our local modes for PbTiC>3 using four 
high-symmetry pointso are more localized than the best 
(three point) RW lattice Wannier functions. 

It is clear that localization cannot be the main qual- 
ity criterion for the construction of local modes. If it 
were, then the ZVR scheme, which uses only one high- 
symmetry k-point to construct a (very localized) lattice 
Wannier function, would be the method of choice. In 
fact, the real test for local mode sets is the degree to 
which they reproduce the energetics of the relevant sub- 
space 7Z. That is, in our case, the degree to which the 
dispersion relations of the effective Hamiltonian 

H e g = i? e ff(Ql> Q2, ■ ■ ■ , Qnu) (8) 
match the real dispersion branches associated with 1Z. 
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In Eq. the variables Qi are the amplitudes of the 
local mode variables, so that H e g can be thought of as 
the "projection" of the complete Hamiltonian into the 
relevant subspace 1Z (which is typically considered as 
energetically decoupled from the rest of the configura- 
tion space of the crystal). The explicit form of H c g will 
depend on the detailed structure of the lattice Wannier 
functions. In particular, the number of distinct coupling 
coefficients (representing the interaction of modes at dif- 
ferent sites) which one should take into account in H c g 
is determined by the spatial extent of the local modes. 

We have constructed effective Hamiltonians for the 
model crystal for each of the three local-mode construc- 
tion schemes discussed above (we obtain the coupling be- 
tween Wo s an d s , by calculating the energy associated 
to the crystal when it is distorted by just these modes). 
The original crystal Hamiltonian involved interactions up 
to fourth nearest neighbors for white atoms. Since the 
local modes involve basically displacements of the central 
white atom, we have kept the same range of interaction 
in H c ff , but now referring of course to fourth nearest local 
modes. This amounts to using ten independent coupling 
coefficients; a larger number of parameters would not be 
reasonable in a practical application. 

ZVR-style local modes are very localized and do not 
couple beyond the fourth neighbor shell, so the consid- 
ered -f/eff includes all the existing interactions. This can 
be seen on panel a) of Figure S the dispersion branches 
computed from H c g match the exact ones at the T point. 
However, H c g gives a poor description of the dispersion 
branches away from T, as it should be expected in view 
of the construction procedure. (Incidentally, the inverse 
of Eq. H leads to Bloch modes which are not normalized 
to unity, except at the T point. The standard analysis 
of -f/eff as given would lead to the low-lying dispersion 
branches in the figure. The higher branches are obtained 
by considering the corresponding generalized eigenvalue 
problem.) Panel b) shows that the i? e ff constructed on 
the basis of RW local modes gives a good qualitative over- 
all description of the dispersion, but fails to match the 
exact branches at the T point (as it should, given that this 
point was used in the construction scheme). The reason 
is that the local modes are more extended, and it is nec- 
essary to include couplings tp-iurther shells (at least up 
to seventh nearest neighbors pj for the match to be essen- 
tially perfect. This means that the RW scheme does npi. 
lead to efficient local modes, in the sense stated above.Ea 
The situation gets worse if more accuracy is needed in 
the overall description of the dispersion branches: the lo- 
cal modes turn out to be more extended, and even more 
coupling terms are needed in H c g. 

In contrast, the local modes constructed following our 
heuristic criterion for localization do exhibit good effi- 
ciency (the dispersion branches do not change much when 
couplings to more than fourth nearest neighbors are in- 
cluded) and provide a very good qualitative match of 
the true branches throughout the BZ (Fig. || c)). (It 



should be noted that our construction scheme does not 
involve any high-symmetry points, hence the offset of the 
branches at T and M. We trade an overall good match for 
perfect accuracy at a few points.) Since a few coupling 
terms are enough to take into account the structure of 
the local modes, and the resulting H c g provides a good 
fit to the true branches, our lattice Wannier functions 
are well suited for the local representation of the relevant 
subspace 1Z. Moreover, they can be improved if needed 
by including more k-points in the integration set, with 
only a minor sacrifice in the compactness of the effective 
Hamiltonian. 

We find these general conclusions to remain valid when 
more complicated interaction models are considered. 

The practical application of our method of local mode 
construction to real materials requires the knowledge of 
the normal modes at general points of the Brillouin Zone. 
This information is easily obtained with modern linear- 
response codes without the need for large supercells. In 
the field of phase transitions, the use of this new scheme 
should enable the study of more complicated situations 
than those considered up to now. Competition of insta- 
bilities associated to different regions of the BZ or com- 
plications derived from anti-crossing phenomena are ex- 
amples in which this method is bound to be useful. On 
the other hand, this work might provide an illustration 
of some of its theoretical underpinnings: the physical in- 
terpretation of .the band representation associated to a 
dispersion bancpj or the symmetry-induced continuity of 
phonon spectraEa are two instances of this. 

IV. CONCLUSIONS 

We have presented a straightforward scheme for the 
construction of very localized lattice Wannier functions, 
with explicit consideration of crystal symmetry. The new 
localization procedure enables a systematic improvement 
in the description of the relevant physics (by simply us- 
ing denser sets of special k-points in the BZ integration) 
while still being quite efficient in regard to the number 
of coupling parameters needed in the effective Hamilto- 
nian. Besides, the present method is straightforward to 
implement. 
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APPENDIX 

Let us study in more detail the equivalence between the 
band symmetry emerging from the transformation prop- 
erties of the Bloch modes-, and the band representation 
associated to a (r, r) setJla This can be done by consid- 
ering the action of G on the (r, r) set of local modes and, 
consequently, on the associated Bloch modes w k . s . 

In order to proceed, we need to formulate the trans- 
formation properties of the modes w° s under the action 
of {-R|v} 6 G r . Since we consider {-R|v} acting on the 
modes themselves and not on their components, we de- 
note symmetry operations by the associated operators 
0{R\v}. We have 

0{R\v}tv°„ = J2 D ls(i^})<h (9) 

h=l 

where D T ({^|v}) is the matrix associated to {-R|v} by 
irrep t. Now, we consider the rest of elements in the 
(r, t) set. They can be mathematically defined as 

< >s := 0{E\n}0{Ri\vi}w° a (10) 

where {-E|n} is a lattice translation and {i?i|vi} is one 
of the d r elements in G/G r , which are chosen so that all 
the Tj := {i?i|vi}r lie in the same cell. The action of any 
{i?|v} e G on an arbitrary 10" s can be decomposed in: 
a lattice translation, a change of the center and a local 
transformation. Mathematically, this is expressed as 

0{R\v}(0{E\n}0{Ri\vi}w? !S ) = 
0{E\{R\v}(n + Ti ) - Tj } OiRjlvj} 0{R\4} w° a (11) 

where {Rj\vj} G G/G r and {-R|v| £ G r are univocally 
determined. Together with Eq. |9j, this expression de- 
fines the band representation and, by using the inverse 
of Eq. ||, it can be written in the basis of Bloch modes. 
We obtain 

0{R\v}u* itS = 

d T 

exp(^ J Rk({ J R|v}r i - Tj )) £ Dj ls ({R\v}) u™ h (12) 

h=l 

By examining the representations this equation defines 
in the high symmetry k-stars, it can be easily checked 
whether the band representation (r, r) is, equivalent to 
the band symmetry we want to describe. c3 

Once a convenient (r, r) set is chosen, Eq. [jj] fixes the 
requirements on Bloch modes so that they lead to sym- 
metry adapted local modes to" s . For pure translations, 
Eq. [l2| reduces to Bloch theorem. Point symmetry de- 
termines the transformation properties of the d T d r Bloch 
modes in each k-point and establishes the relationship 
of these with those in the rest of the k-star. However, 
Eq. [l^ does not determine the form of the M k matrices 
completely. For instance, in a general k-star (G k = {E}) 



no condition is imposed on the choice of Bloch modes in 
a representative k, though, once this is done, the modes 
in the rest of the star are fixed. This is the freedom we 
use in our construction procedure to get the localization 
of the modes. 
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FIG. 1. a) Two-dimensional model crystal, b) Dispersion 
branches for the model crystal assuming a simple harmonic 
model. 

FIG. 2. Lattice Wannier functions for the optical band of 
the model crystal. Only the function transforming as the first 
component of the two-dimensional representation E of 4mm 
is shown on each panel, a) ZVR method using information at 
T; b) Rabe-Waghmare method using information at (P, X, M); 
c) This work, using a three k-point special set. For c), dis- 
placements can be calculated for up to six shells around the 
o atom. 

FIG. 3. Optical dispersion branches for the model crystal 
in the T-M direction. Exact results (lines) are compared to 
those from effective Hamiltonians (squares) corresponding to: 
a) ZVR local modes; b) RW lattice Wannier functions; c) 
This work's local modes (obtained with a 3 k-point set). In 
a) we also present (as circles) the branches obtained if the 
generalized eigenvalue problem is not considered (see text). 



TABLE I. Contribution to the local mode's norm from 
each of the first four shells around the o atom. For RW, 
the complete mode is contained in the first four shells, but it 
is not normalized to unity. For local modes calculated with 
our scheme using 3 and 10 k-points, less than one per cent of 
the norm is outside these four shells. 



Shell 


This work (3 k) 


This work (10 k) 


RW 


1st 


0.8828 


0.8817 


0.8241 


2nd 


0.1100 


0.1099 


0.1617 


3rd 


0.0001 


0.0006 


0.0029 


4th 


0.0009 


0.0014 


0.0021 


Total (1st to 4th) 


0.9938 


0.9936 


0.9908 


Complete mode 


1.0000 


1.0000 


0.9908 
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Fig. la) — Iniguez, Garcia, and Perez-Mato 
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Fig. 2 — Iniguez, Garcia, and Perez-Mato 
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